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The  interaction  between  charge,  heat  and  mass  transfer  occurring  in  SOFCs  is  investigated  applying  a  finite- 
volume-based  SOFC  model.  The  strong  interactions  are  the  consequence  of  the  high  degree  of  integration  of 
different  processes  (chemical/electrochemical  reactions,  diffusion,  heat  and  mass  transfer)  within  SOFCs. 
The  understanding  of  these  interactions  is  a  key  for  the  future  development  and  application  of  SOFCs.  The 
investigation  was  conducted  by  means  of  a  sensitivity  analysis  for  two  different  fuel  gases,  where  one  gas 
features  a  considerable  amount  of  methane  inducing  steam  reforming  reactions  as  additional  disturbance 
factor  in  the  energy  and  mass  balance  system  of  SOFCs.  In  order  to  isolate  the  impact  of  the  varied  model 
parameters  and  the  according  changes  in  the  interactions  of  charge,  mass  and  heat  transfer  from  side 
effects,  the  sensitivity  analysis  was  conducted  at  constant  fuel  utilization.  It  was  found  that  the  impact  of 
different  fuel  gases  on  the  operational  conditions  of  SOFCs  dominates  geometrical  and  material-induced 
phenomena.  The  power  output  was  most  affected  by  the  fuel,  followed  by  the  values  for  the  activation 
polarization  activation  energy  that  reflects  the  employed  electrode  catalysts  activity. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

SOFCs  are  expected  to  play  an  important  role  in  future  power 
generation  due  to  their  high  efficiency.  SOFCs  are  character¬ 
ized  by  a  high  degree  of  integration.  They  can  be  considered  as 
multi-functional  reactors  combing  heat  exchange,  (electro-)  chem¬ 
ical  reactions  (oxidation,  steam  reforming,  and  water  gas  shift 
reaction)  and  ionic  conduction  through  a  membrane.  All  these 
functions  are  strongly  interrelated,  making  SOFCs  highly  complex 
devices.  Assessing  the  interactions  between  heat  and  mass  transfer 
with  chemical  and  electrochemical  reactions  has  become  increas¬ 
ingly  important  for  an  effective  development  of  SOFCs  and  since 
computational  power  improved,  modeling  has  become  the  most 
important  tool  in  this  respect. 

However,  the  complex  mechanisms  of  SOFC  electrochemistry 
are  not  yet  completely  understood  and  are  therefore  mostly  mod¬ 
eled  using  semi-empirical  correlations.  This  introduces  several 
constants  to  the  model  equation  system,  which  have  to  be  deter¬ 
mined  experimentally  in  order  to  allow  for  a  reliable  thermal 
behavior  and  performance  prediction  of  the  SOFC  in  question.  The 
relative  importance  of  model  parameters,  reforming  kinetics  and 
heat  transport  parameters  is  assessed  with  the  presented  sensitiv¬ 


*  Corresponding  author.  Tel.:  +41  56  310  27  06;  fax:  +41  56  310  21  99. 
E-mail  address:  tilman.schildhauer@psi.ch  (T.J.  Schildhauer). 

0378-7753 /$  -  see  front  matter  ©  2008  Elsevier  B.V.  All  rights  reserved. 
doi:10.1016/j.jpowsour.2008.05.044 


ity  analysis  aiming  at  the  impact  on  the  current  density  distribution 
and  local  temperatures  resulting  in  the  predicted  power  output.  The 
analysis  is  carried  out  considering  humidified  hydrogen  and  pre¬ 
reformed  natural  gas.  This  way  the  interrelations  between  charge, 
heat  and  mass  transfer  induced  by  electrochemical  reactions  can  be 
investigated  either  isolated  from  or  including  the  impact  of  steam 
reforming  and  shift  reactions. 


2.  Model  definition 


The  generalized  finite-volume-based  model  was  discussed  in 
depth  in  [1].  Therefore  only  the  general  approach  and  the  model 
extensions,  which  were  conducted,  to  allow  for  this  sensitivity  anal¬ 
ysis  are  discussed  in  this  work. 

The  model  consists  of  an  electrochemical  performance  model 
coupled  to  a  mass  balance  and  an  energy  balance  model.  For  a 
given  fuel  composition  and  operational  voltage,  the  electrochem¬ 
ical  performance  model  calculates  the  current  density  originating 
from  hydrogen  and  carbon  monoxide  conversion  according  to  Eqs. 


(l)-(3). 

H2+02"4>  H20  +  2e~ 

(1) 

CO  +  02~  +*■  C02  +2e~ 

(2) 

02  +4e-  +>  202~ 


(3) 
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Nomenclature 

A  area[m2] 

d  diameter  [m] 

D  diffusion  coefficient  [m2  s-1  ] 

E  voltage  [V] 

fact  activation  energy  for  exchange  current  density  cal¬ 
culation  [J  mol-1  ] 

F  Faraday  constant  [C  mol-1  ] 

I  current  density  [A  nrr 2  ] 

/0  exchange  current  density  [A  nrr2  ] 

I<p  equilibrium  constant 

l  length  [m] 

n  number  of  exchanged  electrons  per  electrochemical 

reaction 

ndot  molar  flow  [mol  s-1  ] 

Nu  Nusselt  number 

p  total  or  partial  pressure  [N  nr2 ] 

Pel  electrical  power  (DC)  [W] 

Qdot  heat  flux  [W  m-1  ] 

rj  reaction  rate  of  reaction  j  [mol  s-1  m-1  ] 

R  resistance  [Q] 

R  ideal  gas  constant  [J  mol-1  K-1  ] 

TaI<  anode  gas  temperature  [K] 

TSI<  solid  structure  temperature  [K] 

Greek  letters 

aan  convective  heat  exchange  coefficient  at  anode 

[W  nr2  K-1  ] 

aca  convective  heat  exchange  coefficient  at  cathode 
[Wm-2  K-1  ] 

P  transfer  coefficient 

y  pre-exponential  factor  for  exchange  current  density 

calculation  [Am-2] 

8  thickness  of  component  [m] 

AH  heat  of  reaction  [J  mol-1  ] 

£  porosity  of  porous  media 

rjact  activation  polarization  voltage  loss  [V] 

77diff  diffusion  polarization  voltage  loss  [V] 

Pohm  ohmic  voltage  loss  [V] 

As  solid  structure  heat  conductivity  coefficient 

[Wm-1  K-1  ] 

Aan  thermal  conductivity  of  anode  gas  [W  m-1  K-1  ] 

Aca  thermal  conductivity  of  cathode  gas  [W  m-1  K-1  ] 

Vy  stoichiometric  coefficient  of  specie  i  in  reaction  j 
r  tortuosity  of  porous  media 

Subscripts  and  superscripts 

an  anode  gas  channel,  anode  electrode 

ca  cathode  gas  channel,  cathode  electrode 

circ  circumferential 

conv  convective 

cross  cross-sectional 

ed  educts 

eff  effective 

equiv  equivalent 

hyd  hydraulic 

op  operational 

prod  products 

SH  sensible  heat 

STR  steam  reforming 

tot  total 

TPB  triple  phase  boundary 


The  voltage  balance  which  has  to  be  satisfied  by  the  calculated 
current  density  values  is  given  in  the  following  equation: 

Pop  =  pNernst,/  —  Pact ,t  ~  Pdiff ,i  ~  Pact,02  —  Pdiff,02  —  Pohm 

with  i  =  H2 ,  CO  (4) 


In  Eq.  (4),  pact  denotes  the  activation  polarization  losses,  which 
are  computed  via  the  Butler-Volmer  equation,  Eq.  (5). 


I  =  Io 


exp 


fBnFp  act 

V  RTSK 


) 


-  exp 


(-(1-/8) 


nFr 7act 
RTSK 


(5) 


The  implicit  Butler-Volmer  equation  involves  the  exchange 
current  density,  J0,  which  is  calculated  for  the  electrochemical  reac¬ 
tions  occurring  at  anode  and  the  cathode  according  to  Eqs.  (6)  and 
(7),  respectively. 


with  i  =  H2,  CO  and  j  =  H20,  CO2 


(6) 


— (Pf)“5“p(^r) 

The  ohmic  voltage  loss  is  usually  the  dominating  effect  in  SOFCs. 
The  calculation  is  straightforward  following  Ohm’s  law,  Eq.  (8). 


Pohm  —  Pequiv^tot  (8) 

Besides,  the  present  model  considers  diffusion  polarization 
losses.  These  account  for  the  different  partial  pressures  at  the 
electrochemical  reaction  sites  compared  to  the  bulk  gas  phase. 
Following  Fields  law,  the  partial  pressures  of  the  reactants  at  the 
electrochemical  reactions  can  be  computed  according  to  the  fol¬ 
lowing  equation: 


P™  =  Pi+ 


{  RTaKI8anran\ 
{  nFDeff£  an  J 


with  -  for  i  =  H2,  CO  and  +  for  i  =  EI20,  C02  (9) 


Knowing  the  bulk  gas  and  the  reaction  site  partial  pressure 
values,  the  diffusion  polarization  losses  are  calculated  via  the  fol¬ 
lowing  equation: 

/'RTSK\  ,  fPiPjFB\  .  L  .  „  m  ,  . 

Pm,i=  \~rf-)  ln  (  — W  )  Wlth  ,=h2-  co  and  J=H20,  C02 

7  '  (10) 


The  mass  balance  model  requires  the  current  density  as  input 
for  the  calculation  of  the  conversion  rates  of  the  electrochemi- 
cally  active  species,  which  are  hydrogen,  water,  carbon  monoxide 
and  carbon  dioxide.  Further,  the  mass  balance  model  includes  the 
calculation  of  reaction  rates  for  all  considered  homogenous  and 
heterogeneous  reactions,  Eqs.  (11 )— (13),  based  on  applied  kinetic 
models  or  on  the  equilibrium  approach. 

CO  +  H20  C02  +  H2  (11) 

C02  +  H2*>  CO  +  H20  (12) 

CH4  +  H20  CO  +  3H2  (13) 

A  literature  review  revealed  contradictory  results  with  respect 
to  applied  kinetic  models  for  the  steam  reforming  reaction  (STR), 
Eq.  (13),  over  nickel-cermet  materials.  The  highest  discrepancies 
were  found  in  the  reaction  order  of  water.  Achenbach  and  Riensche 
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Table  1 

Coefficients  of  the  equilibrium  constant  fit  correlation  for  the  methane  steam 
reforming  reaction 


Coefficient 

Methane  steam  reforming  reaction 

yo 

20.45485 

-294.50512 

ti 

69.62891 

-144.02845 

t2 

182.27607 

a3 

-70.17964 

t3 

693.27912 

found  a  reaction  order  of  zero  [2],  Ahmed  and  Foger  as  well  as  Lee  et 
al.  found  negative  reaction  orders  [3,4]  and  Leinfelder  found  a  pos¬ 
itive  reaction  order  [5].  Drescher  has  shown,  that  all  these  findings 
are  correct  and  a  result  of  the  chosen  operating  conditions  of  the 
experiments,  more  precisely  of  the  chosen  steam-to-carbon  (SC) 
ratio  [6].  Small  SC  yield  positive  reaction  orders,  SC  in  the  order  of  2 
yield  reaction  orders  of  zero  while  high  SC  lead  to  negative  reaction 
orders  of  water.  Drescher  formulated  a  Langmuir-Hinshelwood 
type  kinetic  model,  however  for  temperatures  around  650  °C.  All 
the  mentioned  applied  kinetic  models,  Eqs.  (15)-(18)  [2, 3, 5, 6],  as 
well  as  the  equilibrium  approach,  Eq.  (14),  were  implemented  in 
the  developed  model. 

rCH4STR,eq  =  10000.0  mol  m_1  s-1 

PcoPh2 

kPstrPch4Ph 


1  (bar2 


abs'  VcHaPHjO 


(14) 


fcH4STR,Ach  =  4274.0  mol  m  1  s  1  (bar2bs)  'pch4  exp 

/  -82,000Jmor1  \ 

X  RTj<  ) 

fcH4STR,AhF  =  8542.0  mol  m-1  s”1  (bar2bs)_1pcH45PS°035  exP 
-OS'OOOJmor1^ 


(15) 


-l 


'cH4STR,Lei  =  30.8e  +  10 mol m  s  (bar2bs)  pCh4Ph2o  exp 


-205,  OOOJmol 
RTS  K 


-l 


(17) 


288.52  mol  m_2s_1(bar2.)  Vh4Ph2o 
x  exp(-ll ,  000J  morr/RTsK) 

tCH4STR,Dre  =  .  _  — - i - ~  ~  ~  ~ - i -  (18) 

1  +  16.0barabs_1pcH4  +  0.143  barabs_1pH2o 

x  exp(39,  000Jmor1/RTsR') 

The  equilibrium  approach  is  based  on  the  equilibrium  constant 
the  values  of  which  were  taken  from  [7]  and  fitted  with  Eq.  (19). 
The  fitting  coefficients  are  given  in  Table  1. 

X=y0+A1  exp  (=M)  +*2  exp  (f)  +A3  exp  (f)  (19) 

The  mass  balance  along  the  gas  channels  follows  the  differential 
equation  Eq.  (20). 

d  j 

-£i=yfVjJrj  with  i=H2,  CO,  etc.  and  j=WGS,  STR,  etc.  (20) 


The  energy  balance  model  serves  for  the  calculation  of  the 
effective  temperatures  in  the  solid  structure  and  the  gas  channels. 
It  includes  the  calculation  of  the  purely  convective  heat  transfer 
between  the  gases  and  the  solid  structure  as  well  as  the  conduc¬ 
tive  heat  transport  within  the  solid  structure.  The  convective  heat 
stream  between  the  anode  gas  and  the  solid  structure  is  computed 
following  Eq.  (21). 

Qconv.an  =  aan/drc,an(Pa^  -  Psi^)  with  aan  =  (21) 

“hyd,an 

A  similar  expression  is  formulated  for  the  convective  heat 
exchange  between  the  solid  structure  and  the  cathode  gas,  Eq.  (22). 


Qconv,ca  —  ^ca^circ,ca(^sf^  —  TCK)  with  Oic  a  —  ~j -  (22) 

“hyd,ca 

The  convective  heat  terms  together  with  all  other  heat  source 
and  sink  terms  yield  the  second  order  derivate  differential  equation 
for  the  solid  structure  energy  balance  including  heat  conduction, 
Eq.  (23). 


ASA 


cross 


d2TsK 

dx2 


—  Qconv,ca  —  Qconv,an  —  0sH,ed,an 


(23) 


3.  Sensitivity  analysis 

The  applied  methodology  for  the  sensitivity  analysis  is  that  a 
single  model  parameter  or  a  set  of  interdependent  parameters  is 
varied  at  once  while  all  other  parameters  are  held  constant  at  the 
reference  value.  Unless  otherwise  stated,  the  varied  model  parame¬ 
ter  is  increased  or  decreased  by  25%  and  50%  based  on  the  respective 
reference  value. 

The  operational  voltage  is  the  only  real  control  variable  of  SOFCs 
as  the  fuel  input  flow  is  actually  a  design  value  which  is  defined  by 
the  targeted  fuel  utilization  and  the  total  current.  Further,  the  input 
airflow  follows  from  the  targeted  mean  cell  temperature  and  the 
maximum  allowed  solid  temperatures.  In  order  to  conform  to  this, 
throughout  the  whole  sensitivity  analysis  the  fuel  utilization  and 
total  current  is  held  constant  by  adjusting  the  operational  voltage. 
This  way  the  impact  of  the  varied  model  parameter  can  also  be  bet¬ 
ter  isolated  from  side  effects  being  an  indirect  consequence  of  the 
parameter  variation.  For  instance,  increasing  the  ohmic  resistance 
with  constant  operational  voltage  yields  lower  fuel  utilization.  As  a 
consequence  of  the  reduced  amount  of  electrochemical  reactions, 
the  mean  cell  temperature  drops,  causing  higher  ohmic  losses  and 
so  on.  This  positive  feedback  effect  leads  to  an  overestimation  of  the 
impact  of  the  increased  ohmic  loss.  It  is  not  always  obvious  what 
side  effects  could  possibly  distort  the  impact  of  a  certain  model 
parameter. 

3.1.  Reference  case 

The  values  of  the  varied  model  parameters  are  given  in  Table  2 
for  the  reference  case.  The  model  parameter  used  in  the  prediction 
of  diffusion  losses  were  taken  from  Selimovic  [8],  because  they  rep¬ 
resent  a  porous  material  with  average  properties  considering  the 
technically  possible.  The  temperature-dependent  conductivities 
of  standard  anode,  electrolyte,  cathode  and  ceramic  interconnect 
materials  were  measured  by  Bossel  who  then  derived  the  correla¬ 
tions  proposed  in  [9]  and  used  in  this  work.  The  exchange  current 
density  for  the  activation  polarization  loss  calculation  was  com¬ 
puted  using  Eqs.  (6)-(7)  [10].  The  required  activation  polarization 
model  parameters  were  estimated  by  Costamagna  by  fitting  the 
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Table  2 

Electrochemical  loss  model  input  data  of  BMT 


Model  input  data 

Unit 

Value 

Source 

Activation  polarization 

H2  oxidation  activation  energy 

J  mol-1 

120000.0 

[10] 

H2  oxidation  pre-exponential  factor 

Am-2 

290000000.0 

CO  oxidation  activation  energy 

J  mol-1 

120000.0 

[15] 

CO  oxidation  pre-exponential  factor 

Am-2 

207000000.0 

02  reduction  activation  energy 

J  mol-1 

120000.0 

[10] 

02  reduction  pre-exponential  factor 

Am-2 

700000000.0 

Ohmic  polarization 

Electric  conductivity  of  anode 

ft-1 

T-Dependent 

Ionic  conductivity  of  electrolyte 

ft-1 

T-Dependent 

[9] 

Electric  conductivity  of  cathode 

ft-1 

T-Dependent 

Electric  conductivity  of  ceramic  IC 

ft-1 

T-Dependent 

Diffusion  polarization 

Porosity  of  anode 

- 

0.5 

Tortuosity  of  anode 

- 

3.0 

Average  pore  diameter  of  anode 

m 

1.0E-6 

[8] 

Porosity  of  cathode 

- 

0.5 

Tortuosity  of  cathode 

- 

3.0 

Average  pore  diameter  of  cathode 

m 

1.0E-6 

results  obtained  from  the  model  of  an  electrolyte-supported  pla¬ 
nar  cell  to  experimental  voltage-current  curves  and  temperature 
profiles  [10].  Summing  up,  the  employed  model  parameters  are 
considered  valid  for  the  quantitative  analysis  presented  in  this 
work. 

The  sensitivity  analysis  is  carried  out  for  two  fuel  gases  defined 
in  an  International  Energy  Agency  (IEA)  task  aiming  at  the  definition 
of  a  benchmark  test  (BMT)  for  SOFC  models  (Table  3)  [11  ].  The  IEA  1 
gas  allows  the  explicit  investigation  of  the  impact  of  varying  model 
constants  on  the  electrochemical  performance.  The  IEA  2  gas  is  con¬ 
sidered  to  determine  the  coherences  when  also  internal  reforming 
reactions  take  place.  The  STR  reaction  rates  are  computed  via  Eq. 
(15).  The  cell  geometry  is  planar  and  equal  to  that  of  the  IEA  BMT 
(Table  4).  The  considered  flow  configuration  is  co-flow.  Operational 
conditions  for  the  sensitivity  analysis  are  given  in  Table  5.  The  mate¬ 
rial  property  parameters  used  in  the  energy  balance  calculations  are 
given  in  Table  6. 

Table  7  summarizes  the  integral  model  response  values  for  the 
reference  case  and  both  considered  fuel  gases.  Fig.  1  depicts  the 
predicted  temperature  profiles  of  the  solid  structure  for  the  ref¬ 
erence  case  and  both  considered  fuel  gases.  It  can  be  seen  that 
at  the  cell  inlet  (x-axis:  cell  length  equals  zero)  the  temperature 
predicted  for  the  fuel  gas  IEA  1  is  approximately  65  K  higher  than 
for  the  fuel  gas  IEA  2.  This  can  be  explained  with  the  endother- 
mal  STRs  taking  place  with  the  methane  containing  fuel  gas  IEA 
2.  At  the  cell  outlet  (x-axis:  cell  length  equals  0.1  m),  the  differ¬ 
ence  between  the  two  temperature  distributions  is  decreased  to 
approximately  40  K.  The  difference  between  the  temperature  pro¬ 
files  of  the  two  fuel  gases  is  represented  by  the  curve  “Temperature 
difference”. 


Table  4 

Geometrical  input  data  of  BMT 


3.2.  Variation  of  activation  polarization  parameters 

In  the  model,  the  electrochemical  conversion  of  hydrogen  and 
carbon  monoxide  were  assumed  to  occur  in  parallel.  Each  of  the 
conversion  processes  is  described  with  an  activation  energy  value 
and  a  pre-exponential  factor.  These  model  parameters  are  inter¬ 
dependent  as,  e.g.  increasing  only  the  activation  energy  of  the 
hydrogen  conversion  leads  to  a  compensating  carbon  monoxide 
current  and  vice  versa.  To  avoid  this,  the  activation  energy  values 
of  both  conversion  processes  are  varied  at  once.  The  same  applies 
for  the  pre-exponential  factors. 


Table  5 

Operational  conditions  of  BMT 


Model  input  data 

Unit 

Value 

Source 

Targeted  mean  current  density 

Am-2 

3000.0 

Targeted  fuel  utilization 

% 

85.0 

Air-to-fuel  ratio 

- 

7.0 

Fuel  gas  inlet  temperature 

K 

1173.15 

[11] 

Cathode  gas  inlet  temperature 

I< 

1173.15 

Ambient  temperature 

K 

293.15 

System  pressure 

harabs 

1.01325 

Table  6 

Energy  balance  input  data  of  BMT 


Model  input  data 

Unit 

Value 

Source 

Heat  conductivity  of  anode 

Witt1  K1 

2.0 

Heat  conductivity  of  electrolyte 

2.0 

Heat  conductivity  of  cathode 

2.0 

[11] 

Heat  conductivity  of  interconnect 

2.0 

Nusselt  number 

- 

4.0 

Table  3 

Fuel  and  cathode  gas  compositions  of  the  BMT 
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Table  7 

Model  response  values  for  the  reference  case  for  IEA  1  and  IEA  2  gas 


Response  value 

Unit 

IEA  1 

IEA  2 

Maximum  solid  temperature  gradient 

Knurr1 

2.2 

2.0 

Maximum  solid  temperature 

1327.9 

1286.9 

Minimum  solid  temperature 

K 

1212.0 

1146.3 

Mean  solid  temperature 

1282.7 

1221.9 

Operational  voltage  for  85%  fuel  utilization 

V 

0.683 

0.624 

Cell  lenth  [m] 
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Fig.  1.  Solid  temperature  distributions  and  difference  for  the  reference  case  for  IEA 
1  and  IEA  2  gas. 
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Fig.  2.  Sensitivity  of  current  density  distribution  towards  anode  activation  energy. 


3.2 A.  Activation  energy 

Fig.  2  shows  the  sensitivity  of  the  current  density  distribution 
towards  variation  of  the  anode  activation  energy.  For  the  IEA  1  gas, 
the  maximum  current  is  produced  at  the  cell  inlet.  This  is  due  to  the 
high  hydrogen  partial  pressure  of  the  IEA  1  gas.  From  thereon  the 
current  density  constantly  decreases  along  the  cell  length  because 
of  the  ongoing  fuel  conversion  and  accordingly  decreased  educt 
partial  pressure.  In  contrast,  the  IEA  2  gas  yields  the  maximum 
educt  partial  pressures  in  the  region  where  the  methane  STRs  are 
complete.  Flence,  the  current  density  is  not  constantly  decreasing 
along  the  cell  length,  but  is  somewhat  evenly  distributed  with  a 
wide  and  flat  peak. 

Decreasing  the  anode  activation  energy  means  to  increase  the 
activity  towards  electrochemical  reactions  yielding  higher  current 
density  values  in  the  cell  inlet  region  for  both  gases.  The  increased 
fuel  consumption  in  the  cell  inlet  region  causes  a  considerable  cur¬ 
rent  density  decrease  at  the  cell  outlet  for  the  IEA  1  gas.  However, 
for  the  IEA  2  gas  the  current  density  at  the  cell  outlet  remains  con¬ 
stant,  as  the  fuel  consumption  in  the  cell  inlet  region  is  dominated 
by  the  hydrogen  production  via  STRs.  This  leaves  more  hydrogen 
over  for  the  latter  parts  of  the  cell.  Shortly  spoken,  the  less  hydro¬ 
gen  is  produced  within  the  cell  via  steam  reforming,  the  more  does 
a  decrease  of  the  activation  energy  lead  to  an  increased  current  den¬ 
sity  at  the  cell  inlet  which  is  coupled  to  a  current  density  decrease 
at  the  outlet. 

Table  8  shows  that  for  both  fuel  gases,  the  increase  of  the  anode 
activation  energy  leads  to  an  increase  of  the  temperature  extremes, 
the  maximum  temperature  gradient  and  the  mean  cell  tempera¬ 
ture.  This  is  the  result  of  the  boundary  condition  of  keeping  the 
fuel  utilization  constant,  which  requires  higher  temperatures  due 
to  the  lower  activity  towards  electrochemical  reactions  when  the 
activation  energy  value  is  increased. 

The  operational  voltage  to  meet  the  defined  fuel  utilization 
shows  an  asymmetric  behavior.  It  increases  only  slightly  for  lower 
values  and  strongly  decreases  for  higher  activation  energies  (Fig.  3). 
This  can  be  explained  by  the  investigation  of  the  variation  intervals 
of  the  operational  voltage  for  different  activation  energy  values. 

In  the  reference  case,  the  ohmic  losses  are  dominant;  hence 
the  value  of  activation  losses  is  small.  Decreasing  the  activation 
energy  can  at  the  most  lead  to  a  decrease  of  the  activation  losses 
close  to  zero.  Consequently  the  maximum  increase  of  the  opera¬ 
tional  voltage  is  limited  to  the  activation  loss  value  of  the  reference 
case.  In  contrast,  the  activation  losses  become  rapidly  dominant  by 
increasing  the  activation  energy,  which  is  part  of  an  exponential 
expression.  The  maximum  increase  of  the  activation  losses  is  lim¬ 
ited  by  the  value  of  the  operational  voltage  of  the  reference  case. 
Increasing  the  activation  energy  can  therefore  at  the  most  lead  to 
a  decrease  of  the  operational  voltage  close  to  zero.  The  different 
magnitude  of  the  two  described  variation  intervals  explains  the 


Table  8 

Model  response  values  for  different  anode  activation  energy  values  for  IEA  1  and  IEA  2  gas 


Parameter  value 

Maximum  dT  (Knurr1) 

Cell  outlet  temperature  (I<) 

Cell  inlet  temperature  (I<) 

Mean  cell  temperature  (I<) 

Operational  voltage  (V) 

IEA  1  gas 

-50% 

2.1 

1311.6 

1208.1 

1274.2 

0.750 

-25% 

2.1 

1312.6 

1208.4 

1274.5 

0.746 

Reference 

2.2 

1327.9 

1212.0 

1282.7 

0.683 

+25% 

3.0 

1384.7 

1223.9 

1320.3 

0.446 

+50% 

3.9 

1449.5 

1237.3 

1364.5 

0.175 

IEA  2  gas 

-50% 

1.7 

1264.2 

1142.9 

1207.1 

0.717 

-25% 

1.7 

1265.8 

1143.2 

1208.4 

0.710 

Reference 

2.0 

1286.9 

1146.3 

1221.9 

0.624 

+25% 

2.8 

1345.4 

1150.0 

1255.2 

0.376 

+50% 

3.7 

1408.4 

1153.9 

1290.6 

0.107 
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Fig.  3.  Difference  between  the  predicted  power  output  for  varying  anode  activation  Fig.  5.  Sensitivity  of  current  density  distribution  towards  cathode  activation  energy, 
energy  values. 


observed  asymmetric  behavior  of  the  operational  voltage  for  varied 
activation  energies. 

It  is  important  to  point  out  that  SOFCs  are  usually  not  operated 
at  voltages  below  0.5  V.  The  predicted  voltage  values  for  increased 
activation  energies  indicate  that  in  reality  a  fuel  cell  employing 
such  an  anode  catalyst  cannot  reach  85%  fuel  utilization.  For  the 
lower  activation  energy  values,  the  power  output  increases  by  a 
maximum  of  15%  and  decreases  by  80%  for  the  higher  values. 

Fig.  4  depicts  the  sensitivity  of  the  difference  between  the 
temperature  profiles  of  the  two  fuel  gases  towards  the  anode  acti¬ 
vation  energy.  The  local  temperatures  are  strongly  influenced  by 
the  locally  produced  power,  which  depends  on  the  prevalent  cur¬ 
rent  density  and  the  operational  voltage,  see  also  Table  8.  For  a 
fixed  operational  voltage,  high  current  densities  usually  lead  to 
high  temperatures.  Flowever,  endothermal  reactions  such  as  steam 
reforming  counteract  this,  provided  that  the  reaction  rate  is  con¬ 
siderable. 

It  can  be  seen,  that  higher  values  of  the  activation  energy  lead  to 
an  increase  of  the  difference  between  the  predicted  temperatures 
near  the  cell  inlet  and  a  higher  reforming  peak.  This  can  be  traced 
back  to  the  fact  that  the  temperature  at  the  cell  inlet  for  the  fuel  gas 
IEA  2  is  dominated  by  the  STR  and  does  not  change  significantly  by 
varying  the  activation  energy.  In  contrast,  the  temperature  at  the 
cell  inlet  increases  significantly  with  increasing  activation  energy 
for  the  IEA  1  gas.  The  reason  for  this  is  that  the  operational  voltage 
was  strongly  reduced  in  order  to  meet  the  required  overall  fuel  uti¬ 
lization.  As  a  consequence,  the  amount  of  produced  power  is  low 
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Fig.  4.  Temperature  difference  between  the  predicted  temperature  distributions  for 
IEA  1  and  IEA  2  gas  with  varying  anode  activation  energy  values. 


even  though  the  current  density  is  only  slightly  decreased  with 
increasing  activation  energy.  As  less  electrical  power  is  produced 
from  a  still  considerable  amount  of  converted  fuel,  more  sensi¬ 
ble  heat  is  released  from  the  electrochemical  reactions  resulting 
in  higher  local  temperatures. 

Lowering  the  anode  activation  energy  causes  an  increase  of  the 
temperature  difference  towards  the  cell  end.  Similar  to  the  above 
discussed  increase  of  the  activation  energy,  the  cell  inlet  temper¬ 
ature  of  the  IEA  2  gas  remains  almost  equal  to  the  reference  case. 
Flowever  with  low  activation  energy  values,  the  high  amount  of 
produced  power  plus  the  heat  consumed  by  the  endothermal  steam 
reforming  leads  to  a  lower  temperature  for  the  IEA  2  gas  at  the  cell 
outlet  compared  to  the  reference  case.  For  the  IEA  1  gas  the  missing 
steam  reforming  cooling  is  compensated  to  a  smaller  extent  by  low 
current  densities  at  the  cell  outlet  due  to  fuel  depletion.  Hence,  the 
temperature  difference  between  the  two  gases  is  increased  at  the 
cell  end. 

Fig.  5  depicts  the  sensitivity  of  the  current  density  distribu¬ 
tion  towards  the  cathode  activation  energy.  Note  that  the  cathode 
activation  energy  was  not  decreased  by  more  than  20%  as  the  asso¬ 
ciated  high  electrochemical  activity  at  the  cell  inlet  for  both  fuel 
gases  caused  massive  fuel  depletion  effects  at  the  cell  outlet,  which 
resulted  in  numerical  instability  of  the  employed  solver.  For  the 
reference  case  the  anode  activation  losses  are  slightly  higher  than 
the  cathode  activation  losses.  By  decreasing  the  cathode  activation 
energy  the  anode  activation  polarization  becomes  dominant.  As  a 
consequence  of  the  dependency  of  the  anode  activation  polariza¬ 
tion  on  the  educt  partial  pressures,  the  current  density  increases 
near  the  cell  inlet  or  the  end  of  the  reforming  region  due  to  the 
high  educt  partial  pressures  and  decreases  at  the  cell  outlet  due 
to  the  upstream  fuel  consumption.  This  can  be  observed  for  both 
investigated  fuel  gases. 

In  contrast,  an  increase  of  the  cathode  activation  energy  yields 
a  shift  from  anode  to  cathode  activation  polarization  as  domi¬ 
nant  activation  loss  mechanism.  The  current  density  then  rather 
depends  on  the  oxygen  partial  pressure  than  on  the  educt  partial 
pressures  at  the  anode.  The  oxygen  partial  pressure  does  not  change 
significantly  along  the  cell  length  due  to  the  considerably  smaller 
air  utilization  as  compared  to  the  fuel  utilization.  For  the  IEA  gas  1, 
the  high  current  density  values  at  the  cell  inlet  are  inhibited  by  the 
high  cathode  polarization  and  less  hydrogen  is  consumed.  Conse¬ 
quently,  more  hydrogen  is  left  over  for  conversion  in  the  latter  parts 
of  the  cell.  This  results  in  a  very  even  current  density  distribution 
for  high  cathode  activation  energy  values. 

For  the  IEA  2  gas  an  increase  of  the  cathode  activation  energy 
strongly  reduces  the  already  small  current  density  at  the  cell  inlet. 
Similar  to  the  IEA  1  gas,  more  educts  are  left  for  conversion  in  the 
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Table  9 

Model  response  values  for  different  cathode  activation  energy  values  for  IEA  1  and  IEA  2  gas 


latter  part  of  the  cell.  However,  in  contrast  to  the  IEA  1  gas,  a  current 
density  peak  can  be  observed  at  the  cell  outlet  as  additional  hydro¬ 
gen  is  produced  via  STRs  and  a  still  considerable  oxygen  partial 
pressure  is  present  in  this  region.  Shortly  spoken,  increasing  the 
cathode  activation  energy  yields  a  considerable  shift  of  the  typi¬ 
cal  steam  reforming  induced  current  density  peak  for  hydrocarbon 
containing  fuel  gases  towards  the  cell  end.  The  less  hydrocarbons 
are  present  in  the  fuel  gas,  the  more  even  the  current  density  dis¬ 
tribution  becomes. 

Table  9  shows  that  for  both  fuel  gases,  the  increase  of  the  cathode 
activation  energy  yields  an  increase  of  the  cell  outlet  temperature 
whereas  the  cell  inlet  temperature  remains  almost  constant.  This 
results  in  higher  temperature  gradients.  The  cell  inlet  tempera¬ 
ture  stays  constant  because  the  smaller  amount  of  fuel  conversion 
at  the  cell  inlet  is  compensated  by  more  released  sensible  heat. 
This  is,  because  less  electrical  power  is  produced  due  to  the  lower 
operational  voltage. 

The  impact  of  the  cathode  activation  energy  on  the  required 
operational  voltage  to  maintain  the  defined  fuel  utilization  is 
depicted  in  Fig.  6.  Similar  to  the  analysis  of  the  anode  activation 
energy,  it  is  emphasized  that  operational  voltages  below  0.5  V  are 
not  applied  in  reality.  With  a  maximum  of  approximately  5%  power 
output  increase,  the  decrease  of  the  cathode  activation  energy  has  a 
smaller  impact  as  compared  to  the  decrease  of  the  anode  activation 
energy.  This  can  be  traced  back  to  the  emphasized  dominating  role 
of  the  anode  activation  polarization,  which  is  calculated  based  on 
the  reference  parameters  and  hence  inhibits  a  strong  power  output 
increase. 

Increasing  the  cathode  activation  energy  yields  a  power  out¬ 
put  decrease  of  approximately  75%.  In  this  case,  the  reason  for  the 
smaller  impact  is  that  the  cathode  activation  polarization  is  taking 


Fig.  6.  Difference  between  the  predicted  power  output  for  varying  cathode  activa¬ 
tion  energy  values. 


over  the  dominating  role,  weakening  the  impact  of  fuel  depletion 
towards  the  cell  end.  Thus,  the  power  loss,  due  to  higher  cathode 
activation  losses,  is  counteracted  by  lower  anode  activation  losses 
towards  the  cell  end. 

Fig.  7  shows  the  sensitivity  of  the  difference  between  the 
temperature  profiles  of  the  two  fuel  gases  towards  the  cathode 
activation  energy.  Comparing  Figs.  7  and  4  it  could  be  concluded, 
that  the  cathode  activation  energy  has  a  less  strong  impact  than 
the  anode  activation  energy  values  on  thermal  behavior  of  the  fuel 
cell  for  the  two  investigated  fuel  gases.  However,  the  current  den¬ 
sity  distributions  for  varied  cathode  activation  energy  values  reveal 
another  picture. 

The  strong  reforming  peak,  which  was  discussed  for  higher 
anode  activation  energies,  cannot  be  observed  for  the  increased 
cathode  activation  energy  values.  Similar  to  the  anode  activation 
energy  variation,  the  temperature  at  the  cell  inlet  considering  the 
IEA  2  gas  is  dominated  by  the  occurring  STRs  and  does  hardly 
change  with  varying  cathode  activation  values.  Considering  the  IEA 
1  gas  with  higher  cathode  activation  energy,  the  temperature  at 
the  cell  inlet  is  slightly  lower  then  in  the  reference  case.  The  rea¬ 
son  for  that  is  that  the  produced  current  in  this  region  is  strongly 
reduced  compared  to  the  reference  case  due  to  the  low  oxygen  par¬ 
tial  pressure  at  the  cathode  compared  to  the  initial  hydrogen  partial 
pressure  of  the  IEA  1  gas.  However,  a  strong  temperature  decline 
is  inhibited  by  the  low  operational  voltage,  which  results  in  a  fair 
amount  of  released  sensible  heat,  which  in  return  is  sufficient  to 
maintain  the  inlet  temperature  on  a  level  close  to  the  reference 
case.  Hence,  the  temperature  difference  found  at  the  cell  inlet  for 
increased  cathode  activation  energy  values  is  slightly  lower  than 
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Fig.  7.  Temperature  difference  between  the  predicted  temperature  distributions  for 
IEA  1  and  IEA  2  gas  with  varying  cathode  activation  energy  values. 
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for  the  reference  case,  which  in  return  explains  the  missing  steam 
reforming  peak. 

Towards  the  cell  end,  the  increased  cathode  activation  causes 
higher  current  density  values  due  to  less  fuel  depletion  in  the  front 
parts  of  the  cell.  This  can  be  observed  for  both  fuel  gases.  However, 
the  temperature  increase  towards  the  cell  end  for  the  IEA  2  gas 
is  stronger  due  to  the  already  discussed  current  density  peak  near 
the  cell  end.  This  explains  the  slightly  lower  temperature  difference 
between  the  two  investigated  fuel  gases  at  the  cell  end  for  higher 
cathode  activation  energy  values  as  compared  to  the  reference  case. 

As  discussed  for  the  current  density  distribution,  decreasing 
the  cathode  activation  energy  puts  the  anode  polarization  into  the 
dominating  role.  As  the  anode  polarization  losses  are  computed 
using  the  reference  values,  the  temperature  difference  progres¬ 
sion  remains  almost  unchanged  by  the  lowered  cathode  activation 
energy. 

Concluding  it  can  be  said,  that  the  cathode  and  anode  activation 
energy  values  are  model  parameters,  which  have  to  be  handled 
very  carefully.  The  difference  between  the  local  temperature  dis¬ 
tributions  for  different  fuel  gases,  the  predicted  current  density 
distribution  and  power  output  are  strongly  influenced  by  these 
model  parameters.  In  particular,  the  anode  and  cathode  activa¬ 
tion  energy  values  influence  the  cell  performance  in  a  different 
way.  Thus,  it  is  important  to  know  if  anode  or  cathode  activation 
polarization  is  dominant.  Dominating  anode  activation  enhances 
the  influence  of  the  educt  partial  pressure  distribution  leading  to 
high  fuel  conversion  where  educt  partial  pressures  are  high  and  vice 
versa.  In  contrast,  high  cathode  activation  losses  lead  to  a  relatively 
even  fuel  consumption  along  the  cell  due  to  the  almost  constant 
oxygen  partial  pressures  at  the  cathode  resulting  from  the  usually 
low  air  utilization  resulting  from  usually  excessive  cooling  air  flows. 
Concerning  the  local  temperatures,  high  activation  losses  yield  a 
higher  amount  of  released  sensible  heat  partially  compensating  the 
lower  amount  of  converted  fuel  and  vice  versa.  These  general  rules 
can  only  be  applied  to  a  limited  extent  in  case  internal  reforming 
reactions  occur.  These  dominate  the  local  temperature  in  any  case 
in  the  region  where  they  occur. 

3.2.2.  Pre- exponential  factor 

Within  the  scope  of  variation,  the  impact  of  the  anode  activation 
pre-exponential  factor  on  the  predicted  current  density  distribu¬ 
tion  is  small  compared  to  the  above  discussed  activation  energy. 
Increased  fuel  consumption  at  the  cell  inlet  takes  place  for  higher 
pre-exponential  factor  values.  Analogously  to  decreased  activation 
energy  values,  the  higher  fuel  conversion  at  the  cell  inlet  yields 
lower  current  densities  at  the  cell  outlet  induced  by  fuel  depletion. 


Fig.  8.  Difference  between  the  predicted  power  output  for  varying  anode  pre¬ 
exponential  factor  values. 
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Fig.  9.  Temperature  difference  between  the  predicted  temperature  distributions  for 
IEA  1  and  IEA  2  gas  with  varying  anode  activation  pre-exponential  factor  values. 


However,  the  magnitude  of  changes  is  almost  negligible.  The  sen¬ 
sitivity  of  the  predicted  power  for  constant  fuel  utilization  towards 
the  anode  pre-exponential  factor  is  given  in  Fig.  11.  In  contrast 
to  the  activation  energy,  a  decrease  of  the  pre-exponential  factor 
leads  to  decrease  of  power.  This  is  because,  in  Eqs.  (6)  and  (7),  the 
exchange  current  density  linearly  depends  on  the  pre-exponential 
factor  whereas  it  depends  on  the  multiplicative  inverse  of  the  acti¬ 
vation  energy. 

For  lower  anode  pre-exponential  factor  values,  the  power  output 
decreases  by  a  maximum  of  8%,  while  it  increases  by  a  maximum 
of  4%  for  higher  anode  pre-exponential  factor  values.  Similar  to  the 
variation  of  the  activation  energy  values,  the  IEA  2  gas  shows  a 
more  pronounced  response.  The  reasons  for  that  were  discussed 
above. 

The  minor  changes  of  the  current  density  distribution  appear 
for  both  gases  to  almost  equal  rates.  Further  as  seen  in  Fig.  8,  the 
anode  pre-exponential  factor  has  a  comparably  small  impact  on  the 
operational  voltage  or  predicted  power,  respectively,  for  both  gases. 
Together  this  yields  almost  negligible  changes  of  the  temperature 
difference.  The  reason  for  that  is  that  the  exchange  current  den¬ 
sity  linearly  depends  on  the  pre-exponential  factor  in  contrast  to 
the  exponential  dependency  on  the  activation  energy  value.  Sim¬ 
ilar  results  were  found  in  [12].  Fig.  9  depicts  the  sensitivity  of  the 
temperature  distribution  trend  towards  the  anode  pre-exponential 
factor. 

The  impact  of  the  cathode  activation  pre-exponential  factor  on 
the  predicted  current  density  distribution  is  slightly  more  pro¬ 
nounced  then  for  the  anode  activation  pre-exponential  factor.  In 
contrast,  the  opposite  is  the  case  for  the  impact  on  the  predicted 
power  (Fig.  10).  For  lower  cathode  pre-exponential  factor  values, 
the  power  output  decreases  by  a  maximum  of  5%,  while  it  increases 
by  a  maximum  of  2%  for  higher  anode  pre-exponential  factor  val¬ 
ues.  The  reasons  for  that  are  the  same  as  already  discussed  for  the 
cathode  activation  energy  variation.  Higher  cathode  activation  pre¬ 
exponential  factors  put  the  anode  activation  into  the  dominant  role, 
while  lower  values  lead  to  a  more  homogeneous  current  density 
distribution  which  equalizes  fuel  depletion  losses  at  the  cell  ends. 
Fig.  11  depicts  the  low  sensitivity  of  the  temperature  distribution 
trend  towards  the  cathode  pre-exponential  factor. 

Concluding  it  can  be  said,  that  the  cathode  and  anode  pre¬ 
exponential  factor  values  are  model  parameters  which  have  a 
minor  impact  on  the  difference  between  the  local  temperature  dis¬ 
tributions  for  different  fuel  gases,  the  predicted  current  density 
distribution  and  the  predicted  power  output. 
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Fig.  10.  Difference  between  the  predicted  power  output  for  varying  cathode  pre¬ 
exponential  factor  values. 
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Fig.  12.  Sensitivity  of  current  density  distribution  towards  cathode  diffusion  mate¬ 
rial  parameters. 


3.3.  Variation  of  diffusion  polarization  parameters 

The  calculation  of  the  diffusion  losses  needs  as  input  the  com¬ 
position  of  the  gas  mixture  and  the  properties  of  the  porous  media 
through  which  the  specie  of  interest  diffuses.  The  porous  media 
properties  include  the  pore  diameter,  which  is  important  for  Knud- 
sen  diffusion,  and  the  porosity  as  well  as  the  tortuosity.  Besides  the 
gas  composition  and  the  properties  of  the  porous  media,  diffusion 
losses  depend  on  temperature,  the  diffusion  path  length  and  local 
current  density.  The  latter  is  causing  the  concentration  gradient 
between  the  triple  phase  boundary  (TPB)  where  the  electrochemi¬ 
cal  reactions  take  place  and  the  bulk  gas  phase.  The  diffusion  path 
length  has  a  direct  influence  on  the  concentration  at  the  end  of 
the  diffusion  path.  In  particular,  the  educt  partial  pressure  at  the 
TPB  decreases  with  increasing  diffusion  path  length.  Hence,  when 
talking  about  diffusion  losses,  one  has  always  to  keep  in  mind 
what  kind  of  support  design  is  investigated.  For  instance  stan¬ 
dard  tubular  cathode  supported  cells  feature  a  stronger  cathodic 
diffusion  limitation  than,  e.g.  electrolyte-supported  cells  due  to 
their  approximately  20  times  thicker  cathode.  However,  at  low 
and  intermediate  current  densities,  diffusion  losses  are  usually 
small  or  even  negligible.  Nevertheless,  this  sensitivity  study  has 
been  carried  out  considering  an  electrolyte-supported  planar  cell  in 
order  to  investigate  whether  the  variation  of  the  according  material 
properties  could  lead  to  a  change  of  these  functional  relation¬ 
ships. 
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Fig.  11.  Temperature  difference  between  the  predicted  temperature  distributions 
for  IEA  1  and  IEA  2  gas  with  varying  cathode  activation  pre-exponential  factor  values. 


The  variation  of  the  pore  diameter  was  omitted  as  it  simply 
enhances  the  share  of  Knudsen  diffusion  for  smaller  pore  diameter 
values  or  reduces  it  for  larger  pore  diameters.  Instead,  the  material 
parameters  influencing  both  diffusion  mechanisms,  i.e.  the  poros¬ 
ity  and  tortuosity,  were  varied.  For  the  reference  case  the  quotient 
of  the  interdependent  porosity  and  tortuosity  has  the  value  0.166 
assuming  a  tortuosity  of  3  and  a  porosity  of  0.5.  It  was  decreased  by 
70%  to  0.05,  which  corresponds  to  a  porosity  of  0.25  and  a  tortuosity 
of  5  (typical  values  of  so  called  dense  membranes),  and  increased 
by  110%  to  0.35,  which  corresponds  to  a  porosity  of  0.7  and  a  tor¬ 
tuosity  of  2,  in  order  to  cover  the  full  range  of  technically  possible 
values. 

The  current  density  distribution  did  not  change  at  all  when 
anode  material  parameters  were  varied.  With  increasing  porosity 
and  decreasing  tortuosity,  diffusion  limitations  are  reduced  and  dif¬ 
fusion  losses  become  almost  zero.  This  is  however  hardly  noticeable 
as  diffusion  losses  are  already  two  orders  of  magnitude  smaller 
than  all  other  voltage  losses  in  the  reference  case.  With  dense 
membrane  like  properties  of  the  anode  electrode,  the  diffusion  of 
carbon  monoxide  is  noticeably  limited.  However  the  correspond¬ 
ingly  reduced  carbon  monoxide  conversion  is  compensated  by  an 
increased  hydrogen  conversion,  which  is  possible  because  of  its 
very  high  diffusion  coefficient. 

In  contrast,  Fig.  12  gives  the  sensitivity  of  the  current  density 
distribution  towards  the  cathode  material  parameters.  Increasing 
the  porosity/tortuosity  quotient  does  not  show  any  effect.  This  can 
be  explained  by  the  magnitude  of  the  diffusion  losses  for  the  inves¬ 
tigated  operational  conditions,  which  is  about  two  orders  lower 
than  ohmic,  and  activation  losses.  However,  the  current  density 
progression  computed  with  the  dense  membrane  values  shows  a 
considerable  decrease  of  the  current  density  at  the  cell  inlet  for  the 
IEA  1  gas.  This  indicates  that  in  this  region,  the  local  current  density 
is  close  to  the  diffusion  limited  current  density  of  the  cathode.  In 
these  parts  of  a  cell  the  high  local  current  density  results  from  the 
higher  educt  partial  pressures  as  was  already  repeatedly  discussed. 
For  the  IEA  2  gas,  the  same  phenomenon  can  be  observed,  however 
less  strongly  pronounced. 

For  the  increased  porosity  and  decreased  tortuosity,  the  power 
output  remains  constant.  In  contrast,  assuming  dense  membrane 
properties  for  the  anode  yields  a  decreased  power  output  of  0.5%. 
Decreasing  the  cathode  material  properties  quotient  results  in  a 
2.5%  power  output  decrease.  The  reason  for  this  somewhat  stronger 
response  can  be  found  in  the  smaller  effective  diffusion  coefficient 
of  oxygen  compared  to  hydrogen. 

Fig.  13  shows  the  sensitivity  of  the  temperature  distribution 
trend  prediction  towards  the  cathode  diffusion  material  parameter. 
The  parameter  setting  corresponding  to  a  dense  membrane  leads 
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Fig.  13.  Temperature  difference  between  the  predicted  temperature  distributions  FiS- 14-  Temperature  difference  between  the  predicted  temperature  distributions 
for  IEA  1  and  IEA  2  gas  with  varying  cathode  diffusion  material  parameters.  f°r  IEA  1  anc^ IEA  2  §as  different  Nusselt  numbers. 


to  a  slightly  decreased  cell  inlet  temperature  for  the  IEA  1  gas,  as 
less  fuel  is  electrochemically  converted  than  for  the  reference  case 
at  the  cell  inlet.  For  the  IEA  2  gas  the  temperature  remains  on  the 
reference  level  due  to  only  minor  changes  of  the  prevalent  current 
density  and  the  dominating  STRs.  In  sum,  the  temperature  differ¬ 
ence  decreases  at  the  cell  inlet  for  the  dense  membrane  parameter 
setting. 

Concluding  it  can  be  said,  that  for  the  investigated  electrolyte- 
supported  cell  the  anode  diffusion  material  parameters  have  almost 
no  impact  on  the  model  results.  In  contrast,  the  cathode  diffusion 
material  parameters  can  be  set  in  a  way,  that  diffusion  limitations 
become  relevant  to  some  extent.  It  is  important  to  point  out  that 
this  is  especially  the  case  when  other  support  designs  with  long 
current  paths,  at  the  cathode  or  anode,  are  investigated.  Neverthe¬ 
less,  for  the  investigated  support  design  and  operating  conditions 
the  impact  on  the  predicted  power  output  and  the  thermal  behavior 
can  be  neglected. 

3.4.  Variation  of  heat  exchange  parameters 

Our  model  considers  convective  heat  transport  between  the 
gases  in  the  gas  channels  and  the  solid  structure  as  well  as  conduc¬ 
tive  heat  transport  within  the  solid  structure.  As  can  be  seen  in  Eqs. 
(21)  and  (22),  the  convective  heat  flows  depend  on  the  respective 
heat  transfer  coefficients,  the  heat  exchange  surface  area  and  the 
temperature  difference  between  the  gases  and  the  solid  structure. 
The  heat  transfer  coefficient  is  a  function  of  the  Nusselt  number 
(Nu),  the  thermal  conductivity  of  the  considered  gas  mixture  and 
the  hydraulic  diameter  of  the  gas  channel.  In  some  models  the  Nus¬ 
selt  number  is  considered  variable  [13,14],  however  in  our  model, 
Nu  is  assumed  to  be  constant  with  a  value  of  4.  This  value  lies 
between  the  two  special  cases  of  heat  transfer  with  constant  heat 
flux,  Nu  equal  to  4.36,  and  with  constant  surface  temperature,  Nu 
equal  to  3.66,  and  is  therefore  considered  to  be  a  valid  assumption 
[15].  The  sensitivity  analysis,  considering  both  mentioned  values, 
was  carried  out  in  order  to  assess  the  impact  and  the  accuracy  loss 
of  the  presented  assumption.  Nu  values  higher  than  4.36  and  lower 
than  3.66  were  not  considered.  Lower  Nu  values  are  not  likely  to 
be  found  for  the  typical  operating  conditions  of  SOFCs  and  higher 
Nu  values  can  be  found  for  turbulent  flow  conditions,  which  were 
excluded  from  the  developed  model. 

The  current  density  distributions  of  both  investigated  gases  are 
not  changed  by  varying  the  Nu  values.  As  a  consequence,  the  change 
of  predicted  power  output  is  in  the  magnitude  of  0.1%  and  therefore 
negligible.  Fig.  14  gives  the  sensitivity  of  the  temperature  distri¬ 


bution  trend  prediction  towards  Nu.  The  lower  Nu  value  does  not 
have  a  noticeable  impact  on  the  predicted  temperature  distribu¬ 
tion  difference  for  the  two  investigated  fuel  compositions.  However, 
the  increased  Nu  value  yields  slightly  lower  temperatures  at  the 
cell  inlet  for  the  IEA  1  gas  and  higher  temperatures  for  the  IEA  2 
gas,  yielding  an  overall  smaller  temperature  difference  at  the  cell 
inlet. 

The  observed  opposite  trends  of  both  gases  can  be  explained  by 
the  driving  temperature  difference  between  the  solid  structure  and 
the  inlet  gas  temperature.  For  the  IEA  1  gas,  the  solid  temperature 
is  approximately  27 1<  higher  than  the  inlet  gas  temperature.  There¬ 
fore,  the  inlet  gases  are  heated  up  by  the  solid  structure,  which  in 
return  is  cooled  down.  Increased  heat  transfer  therefore  results  in 
increased  cooling  of  the  cell  inlet  region.  In  contrast,  the  strongly 
endothermal  STRs  occurring  for  the  IEA  2  gas  at  the  cell  inlet  lead 
to  solid  structure  temperatures  approximately  30  K  below  the  inlet 
gas  temperature,  i.e.  the  solid  structure  is  heated  up  by  the  inflow¬ 
ing  gases.  Consequently,  increased  heat  transfer  leads  to  increasing 
solid  structure  temperatures,  which  explains  the  observed  opposite 
trends. 

The  energy  balance  of  the  solid  structure  is  given  by  Eq.  (23). 
From  the  equation  it  can  be  derived  that  temperature  of  the 
solid  structure  depends  on  the  source  terms  originating  from 
chemical  and  electrochemical  reactions,  the  convective  heat  trans¬ 
port  between  gases  and  solid  structure  and  the  conductive  heat 
transport.  The  conductive  heat  transport  is  a  function  of  the 
cross-sectional  area  of  the  solid  structure  through  which  heat  is 
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Fig.  15.  Sensitivity  of  current  density  distribution  towards  the  solid  heat  conduc¬ 
tivity  coefficient. 
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Table  10 

Model  response  values  for  different  solid  heat  conductivity  coefficients  for  IEA  1  and  IEA  2  gas 


transported,  the  temperature  difference  within  the  solid  struc¬ 
ture  and  finally  the  effective  solid’s  heat  conductivity  represented 
by  a  material  parameter.  For  the  reference  case,  the  solid’s  heat 
conductivity  coefficient  (A,s)  was  set  to  2Wm-1  K-1,  which  is  a 
typical  value  for  the  ceramic  components  usually  employed  in 
SOFCs.  The  effect  of  Xs  was  investigated  by  increasing  its  value  to 
20  Wm-1  K-1,  representing  the  value  found  for  high-temperature 
stainless  steel  materials  which  are  used  in  some  state-of-the-art 
SOFCs,  and  decreasing  it  to  0.2  W  m-1  K-1 ,  which  covers  the  case  of 
a  material  with  virtually  no  heat  conduction. 

Fig.  15  presents  the  predicted  current  density  distribution  for  the 
investigated  solid  heat  conductivity  coefficients  for  IEA  1  and  IEA  2 
gas.  Higher  Xs  values  yield  a  strong  increase  of  the  current  density  at 
the  cell  inlet  and  corresponding  decrease  at  the  cell  outlet.  In  con¬ 
trast,  lower  Xs  values  show  a  minor  impact  on  the  predicted  current 
density  distribution,  the  power  output  and  the  thermal  behavior  of 
the  modeled  system  considering  different  fuel  gases.  That  makes 
sense  as  the  reference  case  value  of  Xs  is  already  small.  Lowering  it 
further  reduces  the  already  small  conductive  heat  transfer  to  val¬ 
ues  close  to  zero.  Therefore,  the  further  discussion  focuses  on  the 
impact  of  increased  As  values,  which  induce  stronger  changes  to 
the  energy  balance  of  the  cell. 

Table  10  gives  the  cell  temperatures  and  operational  voltage  val¬ 
ues  for  varying  ks  values.  Comparing  the  cell  inlet  temperature  of 
both  fuel  gases  for  the  reference  case  and  the  increased  solid  heat 
conductivity  it  is  found  that  the  IEA  2  gas  leads  to  a  stronger  increase 
by  approximately  5  K.  In  detail,  the  cell  inlet  temperature  for  the  IEA 
1  gas  increases  by  72  K  versus  77  K  for  the  IEA  2  gas  (Fig.  16).  Simi¬ 
lar  effects  can  be  observed  for  the  cell  outlet  temperatures,  which 
increase,  by  18  K  for  the  IEA  1  gas  and  20  K  for  the  IEA  2  gas.  Over¬ 
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Fig.  16.  Sensitivity  of  temperature  profile  towards  the  solid  heat  conductivity  coef¬ 
ficient. 


all,  these  differences  sum  up  to  a  difference  of  the  increased  mean 
cell  temperatures  of  13  K  for  the  IEA  2  gas  (54  K  increase)  as  com¬ 
pared  to  the  IEA  1  gas  (41  K  increase).  The  difference  of  the  thermal 
effects  for  the  two  investigated  gases  is  ascribed  to  the  endothermal 
STRs  occurring  at  the  cell  inlet  for  the  IEA  2  gas,  which  represent 
a  considerable  heat  sink  and  additional  driving  force  for  solid  heat 
conduction.  Increased  ks  values  allow  more  heat  to  be  transferred 
from  the  rear  parts  of  the  cell  to  the  cell  inlet,  resulting  into  consid¬ 
erably  higher  temperature  compared  to  the  reference  case  for  the 
IEA  2  gas. 

In  contrast,  for  the  IEA  1  gas  no  endothermal  reactions  take  place 
at  the  cell  inlet.  Hence,  the  heat  transfer  from  the  cell  outlet  to  the 
cell  inlet  region  is  only  driven  by  temperature  equalization.  The 
more  pronounced  mean  temperature  increase  due  to  the  less  dis¬ 
tinct  sub-cooling  region  for  the  IEA  2  gas  in  return  elucidates  the 
stronger  impact  of  the  solid  heat  conductivity  coefficient  on  the 
predicted  power  output  which  is  depicted  in  Fig.  17. 

The  lower  Xs  values  lead  to  approximately  1%  power  output 
decrease  while  it  increases  by  a  maximum  of  5.5%  for  higher  A.s 
values.  In  this  respect,  it  is  important  to  point  out  that  all  the 
model  parameters  used  for  the  calculation  of  voltage  losses  were 
held  constant.  Keeping  that  in  mind,  the  resulting  dependence  of 
the  power  output  on  the  solid  heat  conduction  coefficient  is  con¬ 
siderable.  Further  it  gives  an  additional  justification  for  the  strong 
efforts  undertaken  by  SOFC  technology  developers  to  use  metallic 
materials  for  the  interconnector  plates  in  planar  cells. 

Despite  the  inhomogeneous  current  density  distribution  shown 
in  Fig.  15,  higher  Xs  values  lead  to  a  more  homogeneous  tempera¬ 
ture  distribution  along  the  cell  (Fig.  16).  Further,  the  dominating 
influence  of  the  STRs  on  the  cell  inlet  temperatures  is  weak- 


Fig.  17.  Difference  between  the  predicted  power  output  for  varying  solid  heat  con¬ 
ductivity  coefficients. 
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Fig.  18.  Temperature  difference  between  the  predicted  temperature  distributions 
for  IEA  1  and  IEA  2  gas  with  varying  solid  heat  conductivity  values. 


ened  as  the  according  heat  requirements  are  covered  by  the 
enhanced  solid  heat  conduction.  The  consequences  of  this  are 
that  the  predicted  temperature  profiles  for  different  fuel  gases 
converge,  resulting  in  a  flattened  temperature  distribution  with¬ 
out  the  steam  reforming  peak  typical  for  the  reference  case 
(Fig.  18). 

Concluding  it  can  be  said  that  the  Nusselt  number  has  almost  no 
impact  on  the  model  results  which  justifies  the  according  assump¬ 
tion  and  shows  that  the  accuracy  gain  by  holding  Nu  variable  is 
most  probably  small.  The  heat  conductivity  coefficient  strongly 
affects  the  predicted  current  density  distribution  and  local  tem¬ 
perature  distribution.  The  current  density  tends  to  be  less  uniform 
for  higher  heat  conductivity  values  whereas  the  opposite  is  the  case 
for  the  temperature  distribution.  Further,  for  high  solid  heat  con¬ 
duction  coefficients,  the  dominating  role  of  the  STRs  occurring  at 
the  cell  inlet  is  weakened.  Hence,  the  differences  in  local  temper¬ 
atures  between  hydrocarbon  containing  gases  and  pure  syngases 
are  reduced.  A  considerable  impact  regarding  the  predicted  power 
output  was  found  despite  the  fact  that  the  model  parameters  used 
for  the  electrochemical  performance  prediction  were  held  constant 
at  their  reference  value.  However,  the  solid  heat  conduction  is  a 
material  parameter,  which  can  be  reliably  determined.  Therefore  is 
it  not  likely  to  introduce  an  uncontrolled  uncertainty  to  the  model 
results. 
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Fig.  19.  Sensitivity  of  methane,  hydrogen  and  carbon  monoxide  distribution  towards 
STR  kinetics. 
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Fig.  20.  Sensitivity  of  current  density  and  Nernst  voltage  distribution  towards  STR 
kinetics. 

3.5.  Investigation  of  applied  kinetics  models 

In  SOFCs,  hydrocarbons  can  be  converted  to  hydrogen  and  car¬ 
bon  monoxide  via  STR  due  to  the  usually  employed  nickel-based 
anode  material.  The  STR  reaction  rate  affects  the  hydrogen  and  car¬ 
bon  monoxide  partial  pressures  as  well  as  the  related  Nernst  voltage 
distribution  and  the  temperatures  where  the  STR  reactions  take 
place.  Thus,  the  STR  indirectly  affects  the  electrochemical  reactions 
occurring  in  SOFCs.  In  order  to  assess  this  impact,  five  different  STR 
kinetic  models,  Eqs.  (14)-(18),  were  investigated. 

The  reaction  order  of  steam  is  controversially  discussed  in  the 
literature  related  to  STR  kinetic  models.  The  reference  case  of  the 
sensitivity  analysis  employs  the  applied  kinetics  power  law  (PL) 
model  proposed  by  Achenbach  and  Riensche  [2],  who  found  a  reac¬ 
tion  order  equal  to  zero  for  steam.  Ahmed  and  Foger  also  formulated 
a  PL  kinetic  model  for  the  STR  [3].  In  contrast  to  Achenbach,  Ahmed 
found  a  negative  reaction  order  for  steam.  Another  PL  model  was 
developed  by  Leinfelder,  who  derived  a  first-order  dependence 
of  the  STR  on  the  steam  partial  pressure  from  his  experimental 
data  [5].  Drescher  chose  a  Langmuir-Hinshelwood  approach  (LH) 
to  describe  the  STR  of  methane  over  nickel  cermet  electrodes  of 
SOFCs  [6].  LH  approaches  feature  variable  effective  reaction  orders 
depending  on  the  prevalent  partial  pressures  and  temperatures  of 
the  according  species.  The  last  investigated  approach  to  predict  the 
STR  is  based  on  the  equilibrium  constant  of  the  methane  STR. 

Fig.  19  gives  the  axial  distribution  of  methane,  hydrogen  and  car¬ 
bon  monoxide  for  the  different  STR  kinetic  models.  The  Achenbach 
model  predicts  the  end  of  the  STR  region  at  about  a  quarter  of  the 
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Fig.  21.  Sensitivity  of  the  temperature  distribution  for  the  IEA  2  gas  towards  STR 
kinetics. 
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Parameter  value 

Max.  dT (Knurr1) 

Cell  outlet  temperature  (K) 

Cell  inlet  temperature  (K) 

Mean  cell  temperature  (K) 

Operational  voltage  (V) 

IEA  2  gas 

Reference 

2.0 

1286.9 

1146.3 

1221.9 

0.624 

Drescher  LH 

2.2 

1288.4 

1140.7 

1223.2 

0.626 

Ahmed  PL 

2.1 

1286.8 

1143.2 

1221.4 

0.624 

Leinfelder  PL 

3.1 

1293.5 

1141.0 

1229.6 

0.630 

Equilibrium 

3.2 

1293.7 

1141.2 

1230.0 

0.631 

total  cell  length.  The  model  proposed  by  Ahmed  predicts  slightly 
higher  reaction  rates.  Therefore,  STR  reactions  are  complete  after 
a  fifth  of  the  cell  length.  Drescher’s  LH  model  yields  reaction  rates 
which  are  about  half  of  those  obtained  from  the  Achenbach  model. 
The  Leinfelder  model  results  are  almost  congruent  with  the  equi¬ 
librium  approach,  indicating  that  the  nickel-cermet  investigated 
by  Leinfelder  was  highly  active.  Note  that  the  considered  diffusion 
limitation  inhibits  the  completion  of  the  STR  reactions  within  the 
first  control  volume  even  though  equilibrium  is  assumed.  The  faster 
the  STR  occurs,  the  higher  is  the  corresponding  hydrogen  and  car¬ 
bon  monoxide  peak  as  the  electrochemical  consumption  is  slower 
than  the  production  via  the  STR  reactions. 

Fig.  20  shows  the  predicted  current  density  and  Nernst  volt¬ 
age  distribution  for  the  different  STR  kinetic  models.  For  faster  STR 
kinetics,  the  current  density  at  the  cell  inlet  increases,  yielding 
slightly  decreased  current  density  values  at  the  cell  outlet  due  to 
fuel  depletion  effects.  The  reason  for  the  higher  current  density  val¬ 
ues  at  the  cell  inlet  observed  for  faster  STR  kinetics  can  be  found  in 
the  predicted  Nernst  voltage  within  the  STR  region.  Fligher  hydro¬ 
gen  and  carbon  monoxide  partial  pressures  yield  higher  Nernst 
voltages.  Further,  the  Nernst  voltage  increases  with  decreasing  tem¬ 
peratures  due  to  the  endothermal  STR  reactions. 

Fig.  21  and  Table  11  show  that  the  cooling  effect  in  the  STR  region 
is  more  pronounced  for  faster  STR  kinetics  as  the  locally  higher  heat 
requirements  cannot  be  fully  covered  via  heat  conduction  and  con¬ 
vection.  However,  towards  the  cell  outlet,  fast  STR  kinetics  yield 
higher  temperatures  than  the  slow  kinetics.  In  sum,  faster  STR 
kinetics  yields  higher  average  cell  temperatures.  Table  11  further 
gives  the  predicted  maximum  temperature  gradients  and  opera¬ 
tional  voltages  to  maintain  85%  fuel  utilization.  As  a  consequence 
of  the  more  pronounced  cooling  effect  and  the  resulting  higher 
cell  temperatures  towards  the  cell  outlet,  the  maximum  temper¬ 
ature  gradients  increase  with  faster  STR  kinetics.  The  increased 
mean  cell  temperatures  allow  for  higher  operational  voltage  val¬ 
ues  due  to  less  activation  and  ohmic  voltage  losses.  However,  the 


Cell  lenth  [m] 

Drescher  LH  Reference  -  Achenbach  PL - Ahmed  PL  -  -  -  ■  Leinfelder  PL  —  Equilibrium] 


Fig.  22.  Temperature  difference  between  the  predicted  temperature  distributions 
for  IEA  1  and  IEA  2  gas  with  different  STR  kinetic  models. 


gained  power  output  increase  is  in  the  order  of  1%  for  the  equi¬ 
librium  model  as  compared  to  the  reference  case  considering  the 
Achenbach  model. 

Fig.  22  depicts  the  temperature  distribution  trend  prediction 
obtained  for  the  different  STR  kinetics  models.  The  faster  the  STR 
reactions  occur,  the  flatter  and  the  nearer  to  the  cell  inlet  is  the  STR 
peak. 

Concluding  is  can  be  said,  that  the  predicted  temperature  distri¬ 
bution  depends  considerably  on  the  employed  STR  kinetics  model. 
Faster  STR  kinetics  give  rise  to  a  stronger  cooling  effect  of  the  STR 
reactions  which  is  counteracted  by  higher  temperatures  in  the  lat¬ 
ter  parts  of  the  cell  resulting  from  higher  current  density  values 
at  the  cell  inlet.  Concerning  the  predicted  temperature  differences 
for  different  gases,  the  STR  kinetics  model  has  also  a  considerable 
impact.  Especially  the  reforming  peak  observed  for  the  reference 
case  almost  disappears  for  faster  STR  kinetics.  In  contrast,  the  power 
output  is  only  slightly  changed.  The  equilibrium  approach,  which 
yields  almost  instant  methane  conversion,  shows  power  output 
values  increased  by  approximately  1%. 

4.  Conclusions 

A  finite-volume-based  SOFC  model  was  used  for  a  sensitivity 
analysis.  The  analysis  was  conducted  to  gain  insights  in  the  inter¬ 
action  between  internal  charge,  heat  and  mass  transfer  processes 
occurring  in  SOFCs  operated  with  different  fuel  gases. 

In  contrast  to  all  other  varied  parameters,  the  activation  energy 
values,  the  solid  heat  conduction  coefficient  and  the  chosen  kinetic 
model  for  the  steam  reforming  showed  a  considerable  impact 
on  the  predicted  temperature  distributions.  However,  the  impact 
of  the  varied  parameters  on  the  difference  between  the  pre¬ 
dicted  temperature  distributions  for  the  investigated  fuel  gases  was 
always  smaller  than  that  of  the  gases  themselves. 

With  respect  to  the  absolute  power  output  value,  only  the 
activation  energies,  which  reflect  the  activity  of  the  employed 
electro-catalysts,  were  found  to  have  a  strong  influence  within  the 
variation  bandwidth.  All  other  investigated  model  parameters  had 
less  then  10%  variance  to  the  reference  case  power  output.  The 
trends  for  the  power  output  considering  different  fuel  gases  were 
hardly  affected  by  the  varied  parameters. 

It  is  concluded  that  the  impact  of  different  fuel  gases  on 
the  operational  conditions  of  SOFCs  dominates  geometrical-  and 
material-induced  phenomena  characterized  by  the  investigated 
model  parameters. 
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